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ABSTRACT 

The interstellar medium in star-forming galaxies is a multiphase gas in which turbulent support is 
at least as important as thermal pressure. Sustaining this configuration requires continuous radiative 
cooling, such that the overall average cooling rate matches the decay rate of turbulent energy into 
the medium. Here we carry out a set of numerical simulations of a stratified, turbulently stirred, 
radiatively cooled medium, which uncover a fundamental transition at a critical one-dimensional 
turbulent velocity of « 35 km/s. At turbulent velocities below ^ 35 km/s, corresponding to 
temperatures below lO^'^K, the medium is stable, as the time for gas to cool is roughly constant as a 
function of temperature. On the other hand, at turbulent velocities above the critical value, the gas 
is shocked into an unstable regime in which the cooling time increases strongly with temperature, 
meaning that a substantial fraction of the interstellar medium is unable to cool on a turbulent 
dissipation timescale. This naturally leads to runaway heating and ejection of gas from any stratified 
medium with a one-dimensional turbulent velocity above ^ 35 km/s, a result that has implications 
for galaxy evolution at all redshifts. 

Subject headings: galaxies: starburst - ISM: evolution - ISM: structure 



1. INTRODUCTION 

The cosmic history of dark matter is relatively sim- 
ple. Starting as weak perturbations in a Gaussian ran- 
dom field, dark matter structures accrete and merge to 
form ever-larger objects, leading to a hierarchical evo- 
lution that is in excellent agreement with a wide range 
of observations {e.g. Fu et al. 2008; Vikhlinin et al. 2009; 
Komatsu et al. 2011). 

On the other hand, the formation of galaxies within 
these potentials is extremely complex. Although gas ini- 
tially collapses along with the dark matter, it soon feels 
pressure support, cools radiatively, forms molecules, and 
experiences a wide range of intricate physical processes. 
Perhaps the most important of these is the formation of 
stars, which liberate vast quantities of energy that cause 
much of the gas to reverse its history, parting ways from 
the dark matter and traveling outwards into the depths 
of intergalactic space. 

Such large-scale outfiows have been observed in star- 
forming galaxies over a wide range of masses and red- 
shifts {e.g. Heckman, Armus & Miley 1990; Lehnert & 
Heckman 1996; Lehnert, Heckman, & Weaver 1999; Mar- 
tin 1999, 2005; Strickland et a/. 2000, 2004; Pettini et 
a/. 2001; Veilleux, Cecil & Bland-Hawthorn 2005; West- 
moquette et a/. 2007, 2009), and have been shown to 
be of crucial importance in reconciling the distributions 
of galaxies and dark-matter halos {e.g. Scannapieco et 
a/. 2001; Benson et a/. 2003). Yet, at the same time, 
galaxy outfiows are very poorly-understood, as the in- 
terstellar medium (ISM) from which they form can not 
be modeled as a single-temperature gas. Rather it is 
both continuously cooling and being stirred by turbu- 
lence, which in star forming regions is likely to be driven 
primarily by massive stars {e.g. McKee & Ostriker 1977). 
In fact, this simultaneous cooling and turbulent energy 



input is so extreme that it often leads to a supersonic 
medium in which turbulent velocities exceed the thermal 
velocities, and turbulent motions simultaneously support 
the disk and compress a fraction of the gas to drive star 
formation {e.g. Elmegreen & Scalo 2004; Mac Low & 
Klessen 2004). 

The extremely short cooling times within the ISM 
make it impossible to model stellar feedback by simply 
adding thermal energy to the gas, but the range of phys- 
ical scales involved does not allow for the direct simula- 
tion of supernova remnants within a galaxy-scale simu- 
lation. To deal with this problem, galaxy outfiow simu- 
lations have suppressed cooHng with an empirical heat- 
ing function (Mac Low et al. 1989; Mac Low & Ferrara 
1999), imposed an artificial temperature fioor (Tomisaka 
& Bregman 1993; Tenorio-Tagle & Muhoz-Tuhon 1998; 
Suchkov et al. 1994; D'Ercole & Brighenti 1999; Strick- 
land & Stevens 2000; Fujita et a/. 2003; 2004), delayed 
cooling by an arbitrary time delay (Gerritsen & Icke 
1997; Thacker & Couchman 2000), or simply added out- 
flowing gas directly on large scales (Navarro & White 
1993; Mihos & Hernquist 1994; Scannapieco et al. 2001; 
Springel & Hernquist 2003). While outflows arise from 
all of these approximations, the approaches are so simpli- 
fled that several studies have suggested that supernovae 
and stellar winds may not be effective at driving outflows 
at all, but rather the primary driver may instead be ra- 
diation pressure on dust (Thompson et al. 2005) or non- 
thermal pressure from cosmic rays (Socrates et al. 2008). 

On the other hand, Scannapieco & Briiggen (2010), 
showed that by using a subgrid model to track the un- 
resolved velocities and length scales of turbulence driven 
by massive stars, one could produce galaxy-scale outflows 
while maintaining radiative cooling throughout the sim- 
ulation. Furthermore, the material ejected in these simu- 
lations was a multiphase distribution that grew out of the 
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simultaneous turbulent heating and radiative cooling of 
the gas. These results, although built on an approximate 
subgrid approach, illustrated that strong, stellar driven 
turbulence indeed had the potential to produce outflows 
with the properties of observed starbursting galaxies. 

Here we examine this possibility more closely by con- 
ducting a series of direct numerical simulations that 
study the behavior of a stratified, turbulently sup- 
ported medium, representing a small patch of a larger 
galaxy. Previous theoretical work on stratified media 
at these scales can be roughly grouped into two cate- 
gories: (i) two-and three dimensional studies that deposit 
the kinetic energy from a large number of supernovae 
within a single low-density superbubble {e.g. Tomisaka 
& Ikeuchi 1986; Mac Low et al 1989; Tenorio-Tagle et 
al 1990; Tomisaka 1998; Stil et a/. 2009); and (ii) three- 
dimensional studies of the evolution of stratified turbu- 
lence driven by the stochastic injection of thermal and/or 
kinetic energy from individual or small groupings of su- 
pernovae {e.g. Korpi et al. 1999; de Avillez & Breitschw- 
ert 2004, 2005, 2007, 2010; Dib et al. 2006; Joung & Mac 
Low 2006; Cooper et al 2008; Shetty & Ostriker 2008; 
Joung et al 2009; Ostriker & Shetty 2011). 

Here we adopt a third approach, abstracting the prob- 
lem to study the long-term stability of a stratified 
medium that is continuously stirred at a fixed scale, and 
radiatively cooled such that the energy input from turbu- 
lent decay is matched to the overall cooling rate. In this 
way our method is more like the ones adopted in stud- 
ies of compressible turbulence in molecular clouds {e.g. 
Stone et al. 1998; Porter et al. 1999; Padoan & Nordlund 
1999; Cho & Lazarian 2003; Cho et al 2003; Padoan et 
al 2004; Kritsuk et al 2007; Benzi et al 2008; Lemas- 
ter & Stone 2008, 2009; Burkhart 2009), although unlike 
these studies we maintain vertical stratification of the 
medium, do not assume isothermal gas, and include the 
atomic chemical reactions that are important in mod- 
eling the larger scale, higher temperature ISM. While 
this more abstract approach limits the physical processes 
probed by our simulations, it nevertheless has three key 
advantages. 

First, it sidesteps the many uncertain issues surround- 
ing the coupling of massive stars with the surrounding 
medium. In particular, the ISM simulations described 
above not only adopt many approximate techniques to 
deposit the energy from massive stars, but most of them 
also include an additional large source of diffuse heating, 
either taken from a purely numerical model that balances 
cooling in the initial configuration {e.g. de Avillez & Bre- 
itschwerdt 2004) or a model that approximates the pho- 
toelectric gas heating caused by the high-energy electrons 
that are removed from dust grains by stellar ultraviolet 
radiation {e.g. Joung & Mac Low 2006). Our simula- 
tions, on the other hand, allow for cooling at all times 
in all locations, and they do not include any additional 
sources of energy beyond that of the turbulent stirring 
itself. 

Secondly, our approach allows us to study the evolu- 
tion of media with a wide variety of velocity dispersions 
without having to attempt to reproduce them by adjust- 
ing several indirect parameters such as diffuse heating 
rate, supernova rate, and the method used to deposit su- 
pernovae energy in the ISM. Instead, as described below, 
we are able to increase the velocity dispersion simply by 



changing the mean density of the medium, as this in- 
creases the overall average cooling rate, which, in equi- 
librium, is balanced by more vigorous turbulence. In- 
deed, recent simulations that attempt to reproduce the 
turbulent ISM through direct supernova modeling find 
that the one- dimensional turbulent velocity dispersion 
changes only weakly with star formation rate, remain- 
ing < 15 km/s over a wide range of parameters (Dib 
et al 2006; Shetty & Ostriker 2008; Agertz et al 2009; 
Joung et a/. 2009). On the other hand, many galax- 
ies are observed to have velocity dispersions exceeding 
30 km/s, and in some cases even exceeding 100 km/s 
(Forster Schreiber et al 2006; 2009; Wright et al 2007; 
van Starkenburg 2008; Epinat et al 2008, 2009, 2010; 
Cresci 2009; Law et a/. 2009; Genzel et a/. 2008; 2011; 
Green et al 2010; Jones et al 2010; Lemoine-Busserolle 
& Lamareille 2010; Lemoine-Busserolle et al. 2010). 

Finally, and perhaps most importantly, our approach 
lets us separate effects arising from the direct deposi- 
tion of hot gas onto the grid from those caused purely 
by the presence of widespread turbulence. This in turn 
allows us to identify for the first time a global instability 
that occurs in turbulently-supported media driven above 
a critical velocity dispersion. 

The structure of this work is as follows. In §2, 
we describe our numerical simulations of stratified, 
turbulently-stirred, radiatively-cooled media. In §3 we 
present the results of these simulations, focusing on the 
physics of a transition that occurs at a critical turbulent 
velocity dispersion. A discussion is given in §4, which 
also explores the implications of our results for galaxy 
outfiows across cosmic time. 

2. METHODS 

To study the stability of turbulently-supported media 
as a function of velocity dispersion, we carried out a suite 
of simulations using FLASH (version 3.2), a multidi- 
mensional hydrodynamics code (Fryxell et al. 2000) that 
solves the Riemann problem on a Cartesian grid with 
a directionally-split Piecewise-Parabolic Method (PPM) 
(Coleha & Woodward 1984; Colella & Glaz 1985; Fryx- 
eh, Miiller, & Arnett 1989). While FLASH is capable 
of adaptive mesh refinement calculations, in each of our 
simulations, the hydrodynamic equations were solved on 
a fixed grid, in a box of size 2H in the x and y directions, 
spanning —AH to AH in the z direction, where H is the 
gravitational scale height of the gas distribution. In our 
fiducial simulations we adopted a 128 x 128 x 512 grid 
such that H/^x = 64. The boundary conditions were 
taken to be periodic in both the x and and y directions, 
and in the z direction we adopted the FLASH "diode" 
boundary condition, which, like the outflow boundary 
condition, assumes a zero normal gradient for all flow 
variables except pressure, but unlike the outflow bound- 
ary condition, does not allow material to flow back onto 
the grid. 

In the z direction we also assumed a gravitational ac- 
celeration given by 



which contains a smooth transition at z = 0, but is nearly 
constant outside of \z\ = a, where we fixed a = H/2. For 
our initial conditions we adopted a constant temperature 
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of Tinit = 2 X lO^K throughout the simulation and a 
density distribution given by 



p{z,t) = poe 



,2)l/2_ 




(2) 



If it were to cool only moderately to lO^K, this distri- 
bution would remain in hydrostatic equilibrium if go = 
klO^K/ (firripH) where k is the Boltzmann constant, rup 
the mass of the proton and for ionized gas fi = 0.6, and 
we adopted this value of go as our fiducial one. However, 
in cases in which the medium is cooled and stirred vigor- 
ously different values of go may be appropriate, and we 
also explored the impact of varying this value, express- 
ing it in terms of Tg^av such that go = kTgj-g,^ / {/j^iripH) . 
This means that, if it is not slowed down by material 
above it, gas will be able to escape from the midplane to 
the simulation boundaries at ±4i^ if its velocity exceeds 

^esc = V'^goi^H) = 108 km/s V^grav/IO^K. 

Our simulations included seven species: H+, H, H~, 
He, He+, He++, and electrons. The gas was assumed to 
be initially 76% ionized hydrogen and 24% singly ionized 
helium by mass, but this mix was able to change in time 
as we implemented turbulent velocity forcing, cooling, 
and atomic chemistry, as described in more detail below. 



2.1. Velocity Forcing 
The fluid equations solved in our simulations are 

Dp 
Dt 

Dpui dP 



= 0, 



Dt dxi 
DpE dPuj 
Dt dxj 
DpXs 
Dt 



'- P9i + Pf iFcoide 



-- pE cool + pEchem, 



(^2+,2)l/2_ 



-pAsRs- 



(3) 

(4) 
(5) 
(6) 



On the left hand side of these equations p{x,t) is the den- 
sity field, Ui{'K,t) is the velocity field in the i direction, 
P(x, t) is the pressure, ^(x, t) is the total internal energy, 
and Xs(x, t) is the mass fraction of species 5, where t and 

X are time and position variables, and = ^ + 

On the right hand side gi is the acceleration of gravity, 
f(x, t) is a turbulent driving force, Fcoid is the fraction 
of T < 1.5 X 10^ gas, pE cool is the energy loss due to 
radiative losses, pEchem. is the change in internal energy 
due to chemical reactions, and pAgRs is the change in 
the mass fraction of species s due to ionizations and re- 
combinations. 

As our goal is to study the evolution of a cooling ISM 
that is purely heated by the decay of turbulence, no hot 
gas was added directly to the simulations at any time. 
Rather we perturbed the velocity field using the FLASH3 
"Stir" package (Benzi et al. 2008) which we modified as 
described in Pan & Scannapieco (2010). To cleanly dif- 
ferentiate between the impact of turbulence from other 
forms of energy input, we avoided adding PdV energy 
to the system directly through our forcing term, and 
assumed that the fiow driving terms, f(x, t), contained 
only solenoidal modes, i.e., V • f = 0. Note however. Pan 



& Scannapieco (2010) found that the energy dissipation 
timescale in supersonic turbulence is actually indepen- 
dent of the fraction of driving energy contained in com- 
pressible modes. 

We took f to be a Gaussian random vector with an 
exponential temporal correlation. This forcing scheme 
with a finite correlation timescale is different from that in 
studies using an infinitesimal correlation timescale with 
independent forcing at each time step {e.g. Lemaster & 
Stone 2008), or an infinite correlation timescale with a 
fixed driving force {e.g. Kritsuk et al. 2007). However, 
the choice for the value of t/ is not very important be- 
cause it does not affect the energy dissipation timescale 
(Pan & Scannapieco 2010). 

As in Pan & Scannapieco (2010) the behavior of 
the fi term can be summarized as {fi{\<i^t)fj{\<i^t')) = 



Vi{k) {Sij - exp 



{t-t') 



where is the forcing 

correlation timescale, taken to be 20% of H over the 
initial sound speed. Unlike Pan & Scannapieco (2010), 
however, fi was multiplied by two important factors. The 
first factor, Fcoid(^), is the fraction of cold, T < 1.5x lO^K 
gas in our simulation as a function of time. This acted 
to increase the level of stirring as gas cooled, causing a 
feedback loop that results in an equillibrium configura- 
tion with Fcooi roughly fixed. For each simulation we 
varied the magnitude of fi such that this fraction was 
^0.1 over the initial quasi-stable phase discussed below. 
This approach is meant to approximate the situation in 
a real galaxy in which excessive cooling would lead to 
rapid star formation and energy input. Furthermore, the 
second term, exp { — [(^^ + a^)^/^ — a]iJ~^}, causes stir- 
ring to be most vigorous in regions with higher initial 
densities. The inclusion of this term approximates the 
concentration of energy input from stars near their sites 
of formation near the midplane of their host galaxies. 

Also unlike in Pan & Scannapieco (2010), the forc- 
ing wave numbers were chosen to be in the range 10 < 
Hk/2'K < 10.25. To characterize the large scale proper- 
ties of the flow, we deflned a forcing length scale, Lf , from 
the forcing spectrum, Lf = / ^Vf{k)dk/ J Vf{k)dk = 
H/10. This length scale is chosen such that the turbu- 
lent driving occurs on scales substantially smaller than 
the vertical scale height, as would be the case in a real 
galaxy. In fact, as discussed further below, even smaller 
driving scales would be expected in nature, but these 
are not possible in our current simulations because they 
would require adding velocity kicks on unresolved scales. 

As a measure of the overall level of turbulence at 
each timestep, we calculated the volume- weighted aver- 
age one-dimensional velocity dispersion of the material 
as 



Etiu(x,,t)^ 
37V 



^1D,V 



(7) 



where i is a counter over all N cells in the simulation, 
and the factor of 3 accounts for the ratio of the observ- 
able velocity dispersion in a single direction (Jid, from the 
unobservable total velocity dispersion in all three direc- 
tions (J3D. As emphasized in Joung, Mac Low, & Bryan 
(2009) the velocity dispersion is really a scale-dependent 
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quantity 



^1D,V 



(t,r) 



(8) 

where the first integral is over the fuh simulation volume, 
V, and the second integral averages the longitudinal ve- 
locity difference over all point separated by a distance r. 
However we find that at all scales greater than or equal to 
the forcing scale afj^ v(^'^) agrees with ct^d v(^) within 
~ 10% and eq. ([7]) represents a robust value that ac- 
curately characterizes the overall properties of the flow. 
Thus we used a similar method to calculate the mass 
averaged one-dimensional velocity dispersion as 



'1D,M 



it) 



(9) 



and the volume- weighted average ID velocity dispersion 
as a function of height 



3J 



(10) 



where j is now a counter over all J cells at a given height 
z. 

In each of our simulations, we find that the volume- 
weighted average rms turbulent velocity as a function of 
height scales with the density profile as n(z)^/^. This 
can be explained because the height dependence of the 
driving force is chosen to be proportional to the density, 
such that the energy input rate per unit mass goes like 
(X pn{z)'^Lf /aiD,v{z^t)^ where the large eddy turnover, 
^/c^iD,v, enters because it is essentially the timescale 
over which the driving forcing is correlated. Equating the 
energy input rate to the turbulent energy dissipation rate 
y/Lf then gives aioiz^t) (x n ^ , which is biased 
towards the midplane where turbulent driving would be 
the most vigorous in a real galaxy 

Finally, we note that McCourt et al. (2011) and 
Sharma et al. (2011) have identified that the behavior of 
vertically stratified medium simultaneously experiencing 
optically-thin cooling and uniform heating only devel- 
ops a multiphase distribution if the the free-fall timescale 
tff = yH/g is much longer than the overall average cool- 
ing time (where the local cooling time is defined by eq. 
[T3l below). As our simulations are carried out in an equi- 
llibrium configuration in which the average cooling time, 
fcooi is approximately equal to the turbulent dissipation 
timescale, tdiss = ^//c^3D,m, and 



_ H lO^K (jiD^M 
fcooi Lf Tgrav 22 km/s ' 



(11) 



is well above one in all of our runs. Note however, that 
while the average cooling time is always much less than 
tff , the local cooling time can vary significantly from re- 
gion to region due to spatial variations of the flow den- 
sity and temperature. Furthermore, the energy dissipa- 
tion in a turbulent flow exhibits strong spatial fluctua- 
tions, a phenomenon know as inter mittency. This leads 
to nonuniform heating rate by turbulent energy dissipa- 
tion (Pan & Padoan 2009, Pan et al. 2009), which gives 
rise to the spatial variations of the temperature. In fact. 
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Fig. 1. — Top: Radiative cooling rate as a function of tempera- 
ture for solar metallicity gas, as is used in our simulations. Bottom: 
Cooling time of the gas as a function of temperature. The solid 
line shows the cooling time as defined as 1.5nkT/ (pEcool)^ ct^d the 
dotted curve shows the impact of recombinations on the cooling 
time as estimated by (1.5nkT -\- nlS.GeV)/ (pEcool)- 

as we shall see in more detail below, it is the variance in 
energy dissipation and cooling times that lead to a global 
instability. 

2.2. Cooling 

The energy input from turbulent driving in our simula- 
tions, as in a real galaxy, is balanced by radiative cooling. 
This is implemented in the optically-thin limit, assuming 
local thermodynamic equilibrium as 



4ooi = (1 - Y)il - Y/2) 



pHt,z) 



(12) 



where £^cooi is the radiated energy per unit mass, p is 
the density in the cell, rrip is the proton mass, Y is 
the helium mass fraction, fi the mean atomic mass, and 
A(T, Z) is the cooHng rate as a function of temperature 
and metallicity. Here we made use of the tables com- 
piled by Wiersma, Schaye, & Smith (2009) from the code 
CLOUDY (Ferland et al. 1998), making the simplifying 
approximations that the metallicity of the gas is always 
solar and that the abundance ratios of the metals always 
occur in solar proportions. As in Gray & Scannapieco 
(2010), subcycling was implemented within the cooling 
routine itself, such that T and A(T, Z) were recalculated 
every time Ecoo\/E > 0.1. This is equivalent to an inte- 
gral formalism that assumes a constant density over each 
hydrodynamic time step (Thomas & Couchman 1992; 
Scannapieco, Thacker, & Davis 2001). In order to give 
the distribution a chance to attain an initial level of tur- 
bulence, turbulence is stirred as if Fcooi = 0.1 and cooling 
is shut off for the first 0.25 H / {50 km/s) of the simula- 
tion. 

In the upper panel of Figure [H we show the radia- 
tive cooling rates used in our models. Here the peak at 
T ^ 2 X lO^K is dominated by line emission from metal 
ions, whose atomic energy levels are easily excited by col- 
lisions at these temperatures. Above this temperature. 
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many atoms become fully ionized and the effectiveness 
of line cooling decreases, while below this temperature 
collisions become too weak to excite most atomic tran- 
sitions. Finally, the strong drop below about 15,000K 
corresponds to the point at which even the high-energy 
tail of the kinetic energy distribution is unable to excite 
atomic hydrogen transitions. Below this level, cooling by 
dust emission and molecular transitions become impor- 
tant, but these are much less efficient and not considered 
in our simulations, which are focused on the evolution of 
atomic gas. 

In the lower panel of this figure, we plot the radia- 
tive cooling time, defined as the energy per unit volume 
divided by the radiative cooling rate per unit volume 

tcooi = 1.5nkT/(p£;cooi), (13) 

which is inversely proportional to density. From 1.5 x 10^ 
to « 2 X lO^K, A is approximately proportional to T such 
that the cooling time is roughly constant as a function 
of temperature. Above 2 x lO^K, on the other hand A 
drops strongly with temperature, meaning that gas that 
is heated to these temperatures takes orders of magni- 
tude longer to cool. Finally, also shown on this plot is 
an estimate of the cooling time that approximately takes 
into account the additional energy released as gas moves 
from ionized to neutral. This shows that this additional 
energy makes a substantial difference in the evolution of 
the gas below 2 x lO^K, an effect we capture in detail 
by taking account of atomic chemistry throughout our 
simulations. 

2.3. Atomic Chemistry 

To handle ISM ionizations and recombinations, we 
used a modified version of the chemistry solver described 
in Gray & Scannapieco (2010). As in that study, we de- 
fined the molar abundance of the ith species as, Ys = 
Xs/As^ where, as in eq. ([6]) As is the atomic number and 
Xs is the mass fraction of species 5, such that conser- 
vation of mass is given by Xs = 1. The chemical 
evolution of each the 7 species in our simulations was 
cast as a continuity equation in the form 



where Rs is the total reaction rate due to all the binary 
reactions of the form i + j ^ 5 + /, defined as 

Rs = ^Y,YjX,j{s,T) - YsY^Xs^{k,T), (15) 

where Xij{s,T) is the rate at which species i and j react 
to form species 5 at a temperature T. 

All reaction rates were taken from Glover & Abel 
(2008). Because of the intrinsic order of magnitude 
spread in these rates, the resulting equations are 'stiff,' 
meaning that the ratio of the minimum and maximum 
eigenvalue of the Jacobian matrix, Jij = dYi/dYj, is 
large and imaginary. This means that implicit or semi- 
implicit methods are necessary to efficiently follow their 
evolution. To address this problem, we arranged the mo- 
lar fractions of the 6 species (excluding e~) into a vec- 
tor and solved the resulting system of equations using a 
^th Qj^^Qj. Kaps-Rentrop, or Rosenbrock method (Kaps 



& Rentrop 1979) as described in more detail in Gray & 
Scannapieco (2010). 

To ensure the stability of the chemistry routine while 
at the same time allowing the simulation to proceed at 
the hydrodynamic time-step, we developed a method of 
cycling over multiple Kaps-Rentrop time steps within a 
single hydrodynamic time step. To do this we estimated 
an initial chemical time step of each species as 

Ys+0.1Yh+ 

'^cherQjS — ^chem • 7 v-'-^J 

where achem = 0.5 is a constant that controls the desired 
fractional change of the fastest evolving species. Note 
that this estimate is offset by adding a small fraction of 
the ionized hydrogen abundance to deal with conditions 
in which a species is very low in abundance but changing 
rapidly. 

Once calculated, these species are compared to each 
other and the rchem,s associated with the fastest evolving 
species is chosen as the subcycle time step, Atgub- If this 
is longer than the hydrodynamic time step, the hydrody- 
namic time step is used instead and no additional subcy- 
cling is carried out. If subcycling is required, the species 
time step is subtracted from the total hydrodynamic time 
step and the network is updated over the chemical time 
step. At this point we implement the change in internal 
energy due to ionizations and recombinations, which is 
given by 

Atsub^^chem = -Na ^ ^^.AF,, (17) 

where Es and AYg are the binding energy and change 
in the molar fraction of species s over Atgub, and Na = 
6.022 X 10^^. Then the cooling routine is called and given 
the new internal energy the temperature is updated, and 
the species time steps are recalculated and compared to 
the remaining hydrodynamic time step. Given a new 
Atsub, the cycle is repeated until a full hydrodynamic 
time step is completed. Finally, in cases in which the gas 
is > lO^K H and He are always fully ionized, avoiding 
the need for matrix inversions. 

3. RESULTS 

3.1. Scaling Relations and Energy Balance 

The timescale of energy to decay in a turbulent cascade 
is given by tdiss ^ L/V, where L and V are the length 
scale and velocity of the largest turbulent eddies, which 
in our simulations are ^ Lf and cr3D,M = 3^/^criD,M re- 
spectively. This means that if we were to increase the 
physical scale of the entire simulation by a factor X, leav- 
ing everything else the same, the timescale for turbulent 
energy decay would also be increased by a similar factor 
X. 

In a steady state configuration, this energy input is 
balanced by the radiative cooling of the gas, which is 
proportional to the mean gas density, n. This means that 
if we take a system in equilibrium and rescale the domain 
by a factor X while at the same time decreasing the mean 
density by a factor X, the overall timescale will increase 
but the system will otherwise remain the same. For this 
reason, rather than tie our simulations to a fixed length 
scale, we report our results in units of and dynamical 
times tdyn = ^/^3D,M, where d-sD,M is the time averaged 
mass-weighted three-dimensional velocity dispersion. 



6 



^ 2 

< 1.8 

I 1.6 

^ 1.4 

tio 1.2 




6 -- 



0.4 
I 0.3 

tn'^ 0.2 

0.1 

I 0.15 
1 0.1 
0.05 




~i — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 



I I I I I I I I I I I I I I I I 



I I I I I I I I I I I I I I I I 




? I I I I I I I I I I I I I I I I 




J I ' I I I I I I I I I I I I I L 



5 10 15 20 

Fig. 2. — Evolution of global quantities in our fiducial run, S34. 
Top Panel: Mass averaged velocity dispersion (solid line) and vol- 
ume averaged velocity dispersion (dotted line) in units of km/s. 
Second Panel: Mass averaged temperature (sold line) and volume 
averaged temperature (dotted line) in units of K. Third Panel: 
Mass fraction of cold (T < 10^- ^K) gas. Bottom Panel: Mass 
fraction of gas leaving the simulation volume. In all panels time 
is expressed in units of tdyn = 3-'^/^if/a-iD,M where aiD,M = 34 
km/s. 



Furthermore, as the decay rate of the turbulence is 
proportional to the driving scale rather than the scale 
height, and our simulations are unable to achieve val- 
ues of H/Lf as large as encountered in real galaxies, for 
any equilibrium configuration the mean gas densities in 
our simulations will be lower than those in a real galaxy. 
On the other hand, as discussed in more detail below, 
the behavior of the interstellar medium is most likely de- 
termined by the range of temperatures achieved in the 
simulation, which in turn is determined by the mean ve- 
locity dispersion. Thus, while we report Hn in each of 
our simulations, we expect aiD,M to be the key obser- 
vational parameter used to compare them against the 
behavior of real galaxies, and so we name each run ac- 
cording to this value. A summary of aiD,M, Hn, and 
the other important parameters from each of our runs is 
given in Table 1. 

3.2. Fiducial Case 

In Fig. [2] we show the evolution of the global mean 
quantities in our fiducial run, S34. After the short ini- 
tial stirring stage, cooling is turned on and the temper- 
ature begins to drop significantly, leading to Fcoid ^0.1 
and an increase in (Jid.m to ^ 45 km/s. Equilibrium is 
then maintained for roughly three dynamical times, over 
which criD,M and criD,v are roughly equal and constant. 
During this time, the mass and volume averaged temper- 
atures Tm and Ty are similar, and approximately equal 
to the postshock temperature for a v = a3D,M shock 
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where /cb is Boltzmann's constant and fi = 0.6 is the 
molecular weight of the gas. 

This quasi-stable phase is fieeting however, as after ap- 
proximately three dynamical times, the volume averaged 
temperature begins to climb dramatically, increasing by 
over an order of magnitude to 3 x 10^ by t/tdyn ^ 6. 
Furthermore, this rapid proliferation of hot gas is accom- 
panied by a rapid ejection of material, and over 10% of 
the gas mass is driven out from the simulation volume 
from t/tdyn = 3 to 6. 

To explore the cause of this dramatic heating, in Fig. 
Owe present gas phase diagrams for times t/tdyn ^ 1 to 
6. To construct these distributions, we partitioned the 
gas into equally-spaced logarithmic bins in temperature 
and density, and computed the total mass contained in 
each bin. Initially, the majority of the gas is at or be- 
low, T^M5 the postshock temperature corresponding to 
the mass averaged velocity dispersion, and only a small 
fraction of the gas is at somewhat higher temperatures. 
As time progresses, however, a secondary peak develops 
in the PDF, which appears first at ~ 3 x 10^ K and then 
moves to higher temperatures and lower densities as time 
progresses, traveling along a constant pressure trajectory 
with T (X n~^. 

The key to understanding this behavior is to compare 
the local cooling time of the gas to the dynamical time. 
To facilitate this comparison, we first write the dynami- 
cal time in the same units as the cooling times shown in 
Fig. 1, that is 



^1D,M 



1.2 X lO^yrcm" 



Hn 35kms ^ 



Mopc^ 



criD,M 



(18) 



(19) 

which for run S34 gives t^yn^ = 7.2 x lO^yrcm ^. This 
is roughly 10 times the cooling time of lO^K gas, which 
is to be expected as, in equillibrium, gas cooling balances 
turbulent energy input, and the dissipation timescale for 
turbulence is smaller than the dynamical time by a factor 
of H/L = 10. 

However, as it is a stochasitic process, turbulent driv- 
ing does not heat the gas uniformly, but rather deposits 
only a small amount of energy in some regions, while 
heating other regions to several times the mean temper- 
ature, an effect is stronger in supersonic fiows, such as 
those in the ISM (Pan & Padoan 2009, Pan et al. 2009, 
Pan & Scannapieco 2011). As tcooi rises sharply above 
^ 2 X lO^K, the portion of the gas that is heated to above 
this value takes several turbulent dissipation timescales 
to cool. This means that before this overheated gas can 
settle back down to the mean temperature, subsequent 
turbulent energy input will heat it even further, causing 
it to expand to maintain pressure equilibrium with its 
surroundings. Both the increased temperature and lower 
density lead to even longer cooling times, which rapidly 
become orders of magnitude larger than the turbulent 
driving time, and even the dynamical time. The result is 
the rapid, runaway heating of the hot, low-density gas, 
which pushes its way out of the galaxy, dragging a sub- 
stantial portion of the total gas mass with it. In other 
words, the ejection of material from the disk is a direct 
consequence of turbulent support at the level assumed in 
this simulation. 

Note however, that whether this ejected gas is able to 
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TABLE 1 
Run Parameters 



Name 


<^1D,M 

(km/s) 


0'3D,M 

(km/s) 


(Mo pc-2) 


(yr cm~^) 


^grav 

K 


if/ Ax 


^final 
(^dyn) 


-^outflow 


S34 


34 


59 


0.061 


7.5 X lO'* 


1 X 10^ 


64 


20 


0.16 


S20 


20 


35 


0.0061 


1.3 X 10^ 


1 X 10^ 


64 


40 


0.02 


S29 


29 


50 


0.018 


2.6 X 10^ 


1 X 10^ 


64 


20 


0.01 


S61 


61 


106 


0.18 


1.2 X 10^ 


1 X 10^ 


64 


10 


0.95 


S59G 


59 


102 


0.18 


1.3 X 10^ 


3 X 10^ 


64 


10 


0.05 


S35HR 


35 


60 


0.061 


7.3 X 10^ 


1 X 10^ 


96 


12 


0.27 



Note that the column depth at a given velocity dispersion will be significantly less than in a real galaxy due to our inability to adopt a 
turbulent driving scale smaller than Lf /H = 0.1. 
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Fig. 3. — Phase diagram showing logarithmic contours of the mass- weighted PDF of the gas distribution in run S34 at times t/tdyn = 
0.98, 2.0, and 3.0 (top row) and t/tdyn = 4.0, 5.0, and 6.0 (bottom row). In each panel the density is in units of the initial mean density, 
the temperature is in units of K, and the horizontal line corresponds to Td-^ as given by eq. (|18|) with a'iD,M = 34 km/s. All contours are 
labeled by their values relative to the PDF bin with the most mass. 



escape from the host galaxy wih depend on the overah 
gravitational potential. In the case of our fiducial simu- 
lation ^'esc = V^^o^H = 108 km/s which is comparable 
to the vertical thermal motions of lO^K ionized gas, re- 
sulting in a substantial mass of outfiowing gas. Note 
that, assuming a Gaussian distribution of random ve- 
locities, since Vq^^/^id — 3.2, less than 0.2% of the gas 
would be accelerated outward above escape velocity in 
the absence of runaway heating. On the other hand, a 



galaxy with a higher escape velocity, corresponding to ei- 
ther a larger vertical simulation volume or a higher value 
of ^grav7 and would drive material out of the midplane 
without unbinding it from the host galaxy, as studied in 
more detail below. 

Note also that at the same time the high-temperature 
runaway occurs in this simulation, the cold fraction 
shown in Fig. [2] increases appreciably. This can be un- 
derstood in terms of a simultaneous runaway occurring 



8 





1 




4 




4 




4 




4 




4 




4 




2 




2 




2 




2 




2 




2 




X 

\ 

N 





























-2 




-2 




-2 




-2 




-2 




-2 




























-4 




-4 




-4 




-4 




-4 


1 1 Ml 1 1 1 ^ 


-4 





fJVv >N 



2 
4 



1 



t I ' 



.V ,-5'^ -I 





2 




1 2 
x/H 





ii 



2.0 
1.2 

0.3 

-0.5 

-1.3 

-2.2 
-3.0 

0.8 
0.6 
0.5 



0.4 CP 

o 



0.3 



0.0 



x/H 



Fig. 4. — Log of the temperature distribution in K (top), the log of the density distribution in units of the initial mean density (center), 
and the neutral hydrogen fraction (bottom) in a vertical slice through run S34 at times corresponding to those in Fig. [S] From left to right 
^/^dyn = 0.98, 2.0, 3.0, 4.0, 5.0, and 6.0, and lengths are labeled in units of H. 
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Fig. 5. — Evolution of global quantities in a run with a low turbu- 
lent velocity dispersion, S20. Panels and lines are as in Fig. [2] and 
time is expressed in units of tdyn = 3-'^/^if/a'iD,M5 where a'iD,M = 
20 km/s. 



in the cold, dense gas. For this gas, the typical turbulent 
velocity is well above the local sound speed, meaning that 
it is susceptible to shocks with relatively high Mach num- 
bers. Furthermore, because of the short cooling times in 
these regions, these will be strong, radiative shocks that 
lead to substantial density enhancements. Thus in Fig. 
[3l we see a pile up of cold material at 10-100 times the 
mean density, even as the hottest gas moves to lower 
densities and ever-increasing temperatures. This means 
turbulent runaway is a two-way process in which inef- 
ficient cooing in the hottest regions is accompanied by 
compression-driven cooling in the coldest regions. 

The development of this turbulent runaway is further 
illustrated in Fig. HI which shows the temperature, den- 
sity, and neutral hydrogen mass fraction in a vertical 
slice during the same times as the phase diagrams in 
Fig. [3l At t/tdyn = 1, temperature contrasts are modest 
with all the gas well below lO^K throughout the galaxy. 
Likewise, the density distribution is similar to the initial 
hydrostatic configuration, and dense clumps of gas are 
seen only near the midplane, accompanied by regions 
of neutral hydrogen. By t/tdyn = 2, isolated patches 
with temperatures approaching lO^K are visible, and hy- 
drogen near the midplane has been dissociated, but the 
overall temperature and density distributions are largely 
unchanged. 

By t/tdyn = 3, on the other hand, a dramatic change 
has begun. Patches of T > 2 x lO^K are visible near 
the midplane, increasing the pressure to the point that 
the dense clumps of gas are rapidly moved outward, and 
neutral hydrogen is found only at high altitudes. From 
t = 4 onward, heating progresses catastrophically, in- 
creasing the temperature of the gas near the midplane 
up to ~ lO^K, and driving dense clumps of gas to the 
edges of the simulation and beyond. 



3.3. Low-energy Case 

In contrast to our fiducial run. Fig. [5] illustrates the 
evolution of the global quantities in a run in which the 
density, and consequently the cooling and turbulent driv- 
ing rates, are reduced by a factor of ten. In this case, 
the average mass-weighted velocity dispersion is only 
^1D,M = 20 km/s and the corresponding post-shock tem- 
perature is T^j^ = 4.2 X lO^K. After an initial cooling 
phase, the mass averaged temperature settles down to 
approximately this value, maintaining it for the 40 dy- 
namical time we simulated, twice as long as our full fidu- 
cial simulation. Likewise, while some fiuctuations are 
seen in the volume averaged temperature and velocity 
dispersion, which are more sensitive to the low density 
regions near the edge of the simulation, even these settle 
down almost completely after 20 dynamical times. The 
cold fraction remains almost constant at 10% from about 
10 dynamical times onwards, and only 2% of the mass 
leaves the simulation volume in the initial rearrangement 
of the medium. Unlike the fiducial case, the medium has 
settled into a stable equilibrium that can be maintained 
indefinitely. 

Again, this behavior can be directly understood from 
the mass weighed PDF, shown in Fig.[6l Unlike the fidu- 
cial case, the gas phase distributions shown in this plot 
falls entirely below the 2x10^ limit above which tcooi rises 
rapidly. Since the cooling time is roughly constant below 
2 X 10^, the cooling time of overheated gas is comparable 
to the average cooling time, which in turn is comparable 
to the turbulent driving time. This means that the gas 
will have sufficient time to cool back down to the mean 
temperature before further turbulent energy input takes 
place. Furthermore, as turbulent velocities are signifi- 
cantly smaller than in the fiducial run, the shocks in the 
cold regions are much weaker, leading to small density 
enhancements that fail to cause a substantial increase in 

-^cool • 

In addition, the chemical reactions included in our sim- 
ulation further stabilize the mixture of neutral and ion- 
ized gas that occurs in this temperature range. This is 
because when cold gas is heated it dissociates, remov- 
ing internal energy from the system and reducing the 
increase in the temperature. Likewise when hot gas is 
cooled it recombines, adding internal energy to the sys- 
tem that counteracts the cooling. 

Profiles of the temperature, density, and He and He+ 
mass fractions in a vertical slice through this simulation 
are shown in Fig. [71 The temperature fiuctuations are 
minor and dense gas is strongly concentrated near the 
midplane at all times, such that it is extremely convec- 
tively stable. At all altitudes, the gas is a mixture of neu- 
tral and ionized material, and the majority of helium is 
singly-ionized. This singly-ionized helium will recombine 
and release internal energy if it is cooled and dissociate 
and extract internal energy if it is heated. In contrast to 
the run with more vigorous turbulence, the galaxy has 
found a profile that is stabilized from significant changes 
over the lifetime of the system. 

3.4. Intermediate and High-Energy Cases 

Next we consider an intermediate model, S29, in which 
the density is 30% of that in the fiducial run but 3 times 
that in the low-energy run. In this case, whose global 
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Fig. 6. — Phase diagram showing logarithmic contours of the mass- weighted PDF of the gas distribution in run S20 at times t/tdyn = 
10, 20, and 30. As in Figure [S] the density is in units of the ini tial mean density, the temperature is in units of Kelvin, and the horizontal 
line corresponds to T^^, for criD,M= 20 km/s as given by eq. (|18|) . 



quantities are plotted in the left column of Fig. [51 the 
mass and volume averaged velocity dispersions quickly 
settle into a roughly constant configuration, which is 
maintained throughout the simulation. Here aiD,M = 29 
km/s, and while there is some indication of a shift in the 
volume averaged temperature to higher values, this in- 
crease is much less than seen in run S34. Furthermore the 
cold fraction is relatively stable and there is very little 
material lost from the simulation volume. 

An inspection of the mass-weighed PDF, shown in the 
top row of Fig. [9l indicates that while a significant frac- 
tion of the gas is driven to lO^K, where the cooling time 
is comparable to the dynamical time of 2.5 x 10^ yr cm~^, 
this overheated gas falls just short of achieving the tur- 
bulent runaway seen in the fiducial run. All in all, this 
run appears to be at the boundary of the transition be- 
tween stable and unstable turbulent support. Thus it is 
likely that aiD,M ^ 35 km/s represents a critical value 
above which turbulent support will become unstable. 

A very high energy simulation, in which the density 
is increased by a factor of 3 from the fiducial run, con- 
firms that aiD,M ^ 35 km/s is a critical threshold. The 
global evolution from this run, S61, is plotted in the cen- 
ter panel of Fig. [51 which shows an unstable evolution 
even more extreme than in the fiducial run. After only 
one dynamical time, the volume averaged temperature 
begins to exceed the mass averaged temperature, and 
within three dynamical times it moves into the highly 
unstable > lO^K range. At this point, as in the fidu- 
cial run, the combination of long cooing times and the 
rapid energy input from turbulent decay cause the gas 
to quickly be heated to even higher temperatures, in this 
case exceeding lO^K. The result is an extremely rapid 
expansion of the medium, and by t/tdyn = 6 over 90% of 
the gas fiows out of the simulation volume. 

Again, each of these phases is clearly visible in the 
mass-weighted PDF, shown in the bottom row of Fig. 
[H Here we see that at t/tdyn = 2.9 a large fraction of 
the gas has collected at ^ 2 x lO^K, well above the aver- 
age post-shock temperature corresponding to the velocity 
dispersion. By t/tdyn = 4.5 the temperature of the over- 
heated material increases by an order of magnitude, mov- 



ing along a n (X line of constant pressure. Finally, 
by t/tdyn = 6.0 a large fraction of the gas has fiowed out 
of the simulation volume, moving the histogram to the 
left. 

Note that even without turbulent runaway, one would 
have expected this gas to escape from the gravitational 
potential, as the turbulent velocity is high enough that 
the average mass-weighted temperature quickly rises to 
T = 3 X lO^K, while the original gravitational potential 
was established to be in hydrostatic equilibrium with T = 
lO^K. Thus we carried out a second, very high energy 
run, S59G, in which all parameters were the same as in 
run S61, but the gravity was increased by a factor of 3 
such that Tgrav = 3 X lO^K, and Vq^c = 187 km/s. 

In this case, for which the global quantities are plot- 
ted in the rightmost panels of Fig. 8, the evolution of the 
mass averaged and volume averaged velocity dispersions, 
and mass averaged and volume averaged temperatures is 
extremely similar to run S61. In particular, the volume 
averaged temperature again begins to exceed the mass 
averaged temperature within a dynamical time, again 
within three dynamical times it moves into the highly 
unstable > lO^K regime, and again it finally moves off 
to higher temperatures approaching lO^K. However in 
this run, this runaway is not able to drive material off 
of the simulation grid, and the escape fraction after 10 
dynamical times is only 5%. 

Fig. [To] contrasts the temperature evolution in a verti- 
cal slice from this run with the evolution of a similar slice 
from the lower gravity run, S61. Note that outputs at 
slightly different t/tdyn are compared in this figure be- 
cause the dynamical time is dependent on the average 
velocity dispersion, which can not be precisely predicted 
before each run. In both cases, the initially slow build- 
up of overheated gas near the galaxy midplane leads to a 
rapid runaway that pushes dense gas several scale- heights 
outward. Only in the lower gravity case, however, is this 
runaway sufficiently powerful to drive the gas out to dis- 
tances greater than 4i^, while in the high gravity case, 
it reaches heights of only ^ 2H above the disk, where it 
remains until we stopped the simulation at 10 dynamical 
times. This means that while the ejection of gas from 
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Fig. 7. — Log of the temperature (top left panels) and density (top right panels) distributions, and the fraction of neutral He (bottom 
left panels) and He+ (bottom right panels) in a vertical slice from run S20 at times t/tdyn = 10, 20, and 30. 



the disk is a direct consequence of the level of turbu- 
lent support, whether or not this gas will escape from 
the galaxy is dependent on the strength of the gravita- 
tional potential. In a small galaxy with high ^m, the gas 
would escape into the intergalactic medium in a powerful 
galaxy outflow. In a large galaxy with high ^m, on the 
other hand, one would expect the presence of a "galactic 
fountain" type circulation flow in which the gas flows up- 
wards from regions with vigorous turbulence like the one 
were have simulated, but then moves laterally at high 
altitude and falls back onto the disk in less turbulent re- 
gions {e.g. Shapiro & Field 1976; Bregman 1980; Houck 
& Bregman 1990; Walters & Cox 2001), 

3.5. Resolution 

Finally, we assess the impact of resolution on our re- 
sults. As Lf/Ax = 6.4 in the simulations described 
above, a lower-resolution comparison test was not pos- 
sible because the turbulent driving scale would be in- 
sufficiently well resolved. However we carried out a 



higher resolution simulation in which all parameters were 
matched to the fiducial (S34) run, but H/Ax was in- 
creased from 64 to 96. In Figure (TTJ the evolution of the 
global quantities from this run, S35HR, are contrasted 
with those from the fiducial case. Here we see that the 
evolution of both temperature and velocity dispersion is 
very similar in both runs, generating a turbulent run- 
away that commences at almost exactly the same time 
and reaches almost exactly the same final temperature. 

On the other hand, the fraction of the gas flowing out 
of the simulation volume is significantly higher in this 
run. On closer inspection of the gas PDF (which is not 
shown here) this appears to be due to the fact that the 
higher resolution run, which is less susceptible to numer- 
ical diffusion, is able to span a slightly larger range of 
temperatures than the fiducial run. While this differ- 
ence has only a minor effect on the average quantities, 
it does lead to a noticeable increase both in the frac- 
tion of cold, T < 1.5 X lO^K gas, and the fraction of 
T > 2 X 10^ gas that participates in the turbulent run- 
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Fig. 8. — Evolution of global quantities in a run with intermediate turbulent velocity dispersion, S29 (left column), a run with very large 
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away. As the numerical resolution increases, there are 
two effects that are likely to broaden the temperature 
distribution in this manner. The first is that numerical 
diffusion becomes weaker with increasing resolution, re- 
ducing the mixing between hot and cold regions. The 
second is that energy dissipation in turbulent flows is in- 
termittent, meaning that the range of dissipation rates 
broadens at smaller physical scales, causing the heating 
rate to have a broader PDF as resolution increases (Pan 
& Padoan 2009, Pan et al. 2009) 

Together these effects lead to a broader temperature 
PDF in run S35HR, leading to a larger amount of gas at 
extreme temperatures, both hot and cold. The increase 
in hot gas fraction, in particular, is enough to expand 
the overall gas distribution somewhat more than in the 
fiducial run, moving a larger overall fraction of the gas 
off of the grid. Thus we conclude that while the temper- 
ature and density structure of the turbulent runaway are 
well resolved in our simulations, the final vertical distri- 
bution is weakly dependent on resolution and strongly 
dependent on the assumed gravitational profile. 

4. DISCUSSION 

Taken together, our simulations imply that while there 
is more to be learned about the details of gas escape from 
the overall gravitational potential, we can nevertheless 
be confident in our identification of a fundamental tran- 
sition in stratified, turbulently-stirred, radiatively-cooled 
media. The question remains, however, if the transition 
we have found in our simplified models directly corre- 
sponds to a similar transition in the much more complex 
ISM distributions that occur in real galaxies. Based on 
both theory and observations, we argue that there are 
many reasons to believe that it does. From a theoretical 
point of view, one can list a number of ISM processes 
that we do not include in our simulations, yet none of 
them is likely to impede the turbulent runaway we see 



above criD,M = 35 km/s. This runaway arises because 
any medium with sustained turbulence will be continu- 
ously heated by turbulent decay, and this thermal energy 
must be radiated away to maintain equilibrium. This 
means that if turbulence is vigorous enough to move a 
fraction of the gas into a regime in which the cooling 
time increases strongly with temperature, the continued 
deposition of thermal energy will heat this gas dramati- 
cally, driving it out of the midplane, and in many cases 
out of the host galaxy itself. The simplicity of this mech- 
anism suggests that more complete ISM models will only 
refine the connection between this processes and galaxy 
outflows. 

Thus while our nonrotating simulations do not include 
the Coreolis force, this is unlikely to play a key role. 
In particular, its impact can be estimated by calculat- 
ing the Rossby number, the ratio of the typical veloc- 
ity to the shear across the region, which is greater than 
1 in cases in which rotation is negligible. In a simula- 
tion offset from the center of the disk by a distance do, 
this is R ^ cr3D,Tdo/(2^'^o)7 where '^o is the rotational 
velocity at do. For a galaxy with a flat rotation curve 
^0 ~ 9odo where go is the gravitational acceleration as 
in eq ([T]), and for solid body rotation, as in a dwarf 
galaxy or near the center of a large galaxy, Vq — 3^o<^o- 
This means R ~ (criD,T/cs,grav)(<^o/^)"'^^^ for a flat ro- 
tation curve, and R ~ (criD,T/cs,grav)(3c^o/^)"'^^^ for 
solid body rotation, where Cs,grav = (^^grav/z^^p)"^^^ = 
37km/s(Tgrav/10^i^)^/^ This means that R ^ \ \i 
were to consider our lowest velocity simulation to take 
place at the center of a galaxy, and ^ 1 for all other 
simulations and choices of do- Furthermore, the Coreo- 
lis force has no effect for velocities perpendicular to the 
disk, the direction of the gas ejection driven by turbulent 
runaway. Rotation may be a mechanism worth explor- 
ing, but it is one that is unlikely to drastically effect the 





Fig. 9. — Top: Phase diagram showing logarithmic contours of the mass- weighted PDF of the gas distribution in intermediate-energy 
run S29 at times t/tdyn = 3.0, 9.1, and 15, as compared to Taj^ with aiD,M = 29 km/s. Bottom: Phase diagram showing logarithmic 
contours of the mass- weighted PDF of the gas distribution in the very high-energy run S61 at times t/t^yn = 2.9, 4.5, and 6.0, as compared 
to Td-j^ with a'iD,M = 61 km/s. Note that to accommodate the more extreme conditions in this run the maximum temperature in this row 
3 X lO^K and the minimum value of lg(?2/?2o) is —3. 



transition identified here. 

Magnetic fields are also not included in our simula- 
tions, although they are known to exist in galaxies, and 
they can contribute to turbulence via the magentorota- 
tional instability (Piontek & Ostriker 2004; 2005), mod- 
ify the flow dynamics, and affect the decay of turbulent 
velocities into heat. On the other hand, simulations of 
supersonic magnetohydrodynamic turbulence (Stone et 
al. 1998; Mac Low et al. 1998) have shown that the en- 
ergy density of the magnetic fields built up in the MHD 
turbulence are typically ^2 — 3 smaller than the tur- 
bulent energy density, a result that is consistent with 
observations (Ferriere 2001). Furthermore, Lemaster & 
Stone (2009) showed that large changes in magnetic field 
strength only have a minor impact on the dissipation rate 
of turbulence. On the other hand, much like gravity, or- 
dered fields have the potential to influence the flow near 
the plane of the galaxy, even in the presence of large ve- 
locities and high temperatures. Thus, magnetic fields are 
likely to lead to interesting changes in the simulations, 
but unlikely to affect the dynamics strongly enough to 
avoid turbulent runaway at high velocity dispersions. 

While heat conduction has the potential to affect tur- 
bulent heating by moving energy from hot rarified regions 
into the cooler clumps, the addition of this processes to 
our simulations is very unlikely to have any effect on our 



conclusions. Although the precise value of the conduc- 
tivity is relatively uncertain and dependent on magnetic 
field structure {e.g. Cowie & McKee 1977; Malyshkin & 
Kulsrud 2001), it is clear that it is a slow, second order 
process on ISM scales. Instead, mixing between hot and 
cold regions is primarily determined by the cascade of 
structures down to the very small scales at which con- 
duction can operate {e.g. Pan & Scannapieco 2010). In 
fact, if anything, as turbulence is driven at marginally- 
resolved scales in our simulations, the presence of un- 
avoidable numerical diffusion is likely to somewhat over- 
estimate heat conduction and mixing relative to the real 
ISM, as evidenced in the increased mixing in our fiducial 
run, S34, relative to the high resolution run, S35IIR. 

Our simulations also do not include self-gravity effects, 
which can both serve to enhance shear-driven turbulence 
and lead to the formation of Jeans unstable molecular 
clouds. While the swing amplifier instability that occurs 
in self-gravitating rotating disks (Goldreich & Lynden- 
Benn 1965; Julian & Toomre 1966) has been shown to 
be important in low density regions near the edges of 
disk galaxies {e.g. Wada, Meurer, & Norman 2002; Kim 
& Ostriker 2007; Koyama & Ostriker 2009), it is unlikely 
to drive turbulence at the levels considered here. On 
the other hand, molecular clouds can drive more vigor- 
ous turbulence by gravitational scattering (Fukunaga & 
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Fig. 10. — Log of the temperature distribution in run S61 (top 
panels) at times t/tdyn = 2.9, 4.5, and 6.0, as contrasted with those 
from a similar run with higher gravity, S59G (bottom panels) at 
times t/tdyn = 3.0, 4.4, and 6.0. Note that the upper temperature 
in these plots is 3 x lO^K, rather than lO'^K as used in plotting 
lower-energy runs in Figs. U) El and [8] 

Tosa 1989; Gammie, Ostriker, & Jog 1991), and much 
more importantly, by the formation of stars. For any 
complete model of ISM evolution, tracking the forma- 
tion of such clouds is essential, which requires not only 
self-gravity, but a full model of the molecular chemistry 
and cooling. On the other hand, as our goal here is not 
to capture star formation in detail, but rather to explore 
the evolution of stratified media as a function of turbu- 
lent velocity dispersion, this process is only of secondary 
importance. Our aim is not to capture how turbulence 
is driven, but rather study how driven turbulence affects 
the subsequent evolution of the ISM. 

Besides driving turbulence, star formation will also 
have two major effects that are not considered in our sim- 
ulations. First it will lead to diffuse gas heating through 
the photoionization of neutral gas {e.g. Wolfire et al 
2003), and the photoelectric gas heating caused by the 
high-energy electrons that are removed from dust grains 
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Fig. 11. — Evolution of global quantities in run S35HR, a high 
resolution version of the fiducial run. Panels and lines are as in 
Figs.OandO and time is expressed in units of tdyn = 3-'^/'^if/aiD,M 
where o'id,m = 35 km/s. In each panel the thick (colored) lines 
are from run S35HR, while the thin black lines are for comparison 
and are taken from the fiducial run, S34. 

by stellar ultraviolet radiation {e.g. Joung & Mac Low 
2006). To the extent to which such heating varies with 
space and time (Parravano et al. 2003) or contributes to 
gas motions through radiation pressure (Thompson et 
al. 2005) it will contribute to the overall level of turbu- 
lence. To the extent to which such heating is uniform, 
however, it will offset the overall level of cooling, which 
in turn will mean that less turbulent driving is needed 
to maintain equilibrium. In this case, a given turbulent 
velocity dispersion will correspond to a larger overall col- 
umn depth, again making the point that Hn in our simu- 
lations is likely to be much less than measured in similar 
systems in nature, and that our results are much better 
interpreted as a function of criD,M- 

Secondly, supernovae and stellar winds arising from 
star formation will not only drive turbulence, but will 
also deposit hot gas directly into the interstellar medium. 
In fact, the overlapping of such hot regions into a su- 
perbubble is traditionally thought of as the mechanism 
by which outfiows are driven {e.g. Tomisaka & Ikeuchi 
1986; Mac Low et al. 1989; Tenorio-Tagle et al. 1990; Mac 
Low & Ferrara 1999; Strickland & Stevens 2000; Fujita 
et a/. 2003). On the other hand, there is considerable 
observational evidence supporting the idea that galaxy 
outfiows are driven instead by the more widespread in- 
put from the massive stars (Strickland et a/. 2004). For 
example, observed outfiow pressure profiles demonstrate 
that mass and energy are ejected relatively smoothly over 
large regions that match up well with region of intense 
star formation (Heckman et al. 1990; Lehnert & Heckman 
1996; Strickland & Stevens 2000; Strickland et al. 2000). 

This suggests that turbulent runaway may in fact be 
the primary mechanism producing the hot gas observed 
to be pouring out of starbursting galaxies, and that the 
fraction of hot gas directly added by SNe is not only un- 
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certain theoretically, but of secondary importance phys- 
ically. In fact, adding this process to our simulations is 
only likely to obscure the main results. The ISM in galax- 
ies with velocity dispersions below criD,M would contain 
hot gas, but this gas would cool and condense within the 
disk relatively quickly. On the other hand, supernova 
material in criD,M disks would either be mixed into the 
ISM if it were deposited within a dense clump, or escape 
to high latitudes if it were ejected into a low-density, 
superheated region. As seen in our more approximate 
galaxy-scale modeling in Scannapieco & Briiggen (2010), 
this would lead to outflows that were more metal enriched 
than the overall ISM, but this difference would be much 
less than in a "superbubble" outflow picture in which 
the metals from supernovae were injected into a single 
hot cavity. However, as our simulations include energy 
input from stars as random turbulent motions its is not 
possible for us to associate metals from stars with partic- 
ular parcels of gas in our present simulations to quantify 
this effect further. Furthermore, the best way to imple- 
ment supernovae into ISM scale simulations is extremely 
unclear, and to date has lead to turbulent velocity disper- 
sions that vary only weakly with star formation rate and 
increase only up to ^ 15 km/s, (Dib et a/. 2006; Shetty 
& Ostriker 2008; Agertz et a/. 2009; Joung et a/. 2009), 
falling far short of the critical velocity dispersion identi- 
fied here. 

On the other hand, with the advent of modern 
integral-field unit (IFU) spectroscopy, galaxies with tur- 
bulent velocities above 35 km/s have been well ob- 
served at a variety of redshifts, with a variety of meth- 
ods. At cosmological distances, Lemoine-Busserolle et 
al. (2010) and Lemoine-Busserolle & Lamareille (2010) 
conducted Adaptive Optics (AO) IFU spectroscopy with 
the VLT/SINFONI instrument, using rest-frame optical 
nebular emission lines to measure average velocity disper- 
sions of 60 — 70 km/s in four z ^ 3 galaxies and velocity 
dispersions of 35 — 100 km/s in eight z = 1.0 — 1.5 galax- 
ies. At Keck, Jones et al. (2010) conducted AO OSIRIS 
IFU spectroscopic observations of six strongly-lensed z = 
1.7—3.1 galaxies, at resolutions of up to 100 pc, obtaining 
Ha velocity dispersions between (TiD,Ha = 50 — 100 km/s. 
Wright et al. (2007) and Law et al. (2009) used AO and 
the OSIRIS IFU to study a sample of unlensed galax- 
ies between z = 1.5 and 3.1 and found criD,HcK between 
50 — 100 km/s in most cases, and 139 ±3.2 in an extreme 
starbursting region of a single galaxy (Q1623-BX543). 
van Starkenburg et al. (2008) used SINFONI at the VLT 
to observe sl z = 2.03 galaxy with cFiB,Ua ranging be- 
tween 30 — 100 km/s, and Epinat et al. (2009) used this 
instrument to observe nine z = 1.2 — 1.6 galaxies with 
criD,Ha ~ 40 - 100 km/s. Finally, the large SINS sur- 
vey to study the dynamics of high-redshift galaxies also 
used SINFONI IFU spectroscopy, and to date has mea- 
sured over twenty z ^ 2 disk galaxies and star-forming 
clumps with criD,Ha = 30 — 100 km/s (Forster Schreiber 
et al. 2006; 2009;'Cresci 2009; Genzel et al. 2008; 2011). 
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In the nearby universe, the Gassendi Ha survey of SPi- 
rals (GHASP) has used spatially-resolved Fabry-Perot 
spectroscopy to observe a sample of 153 isolated disk 
galaxies with Ha velocity dispersions ranging from 20- 
40 km/s (Epinat et al. 2008; 2010). Selecting a sample of 
higher Ha luminosity disks from the Sloan Digital Sky 
Survey, Green et al. (2010) used the SPIRAL and WiFeS 
IFU spectrographs to measure luminosity-weighted Ha 
velocity dispersions on ?^ 2 kpc scales in a sample of 65 
z ^ 0.1 rapidly star- forming galaxies. These velocity 
dispersions ranged from 25 — 100 km/s, and they were 
strongly correlated with star formation rate and uncor- 
related with stellar mass, suggesting massive stars as the 
primary drivers of vigorous turbulence. 

In observed star forming disk galaxies, the total pres- 
sure in the disk is just sufficient to stabilize all per- 
turbations on scales too small to be stabilized by rota- 
tion, meaning that the Toomre parameter Q = {afj^ + 
c^)V2^/7rGE^ ^1-2, where K is the epicyclic fre- 
quency and Tig is the gas density per unit area (Leroy et 
al. 2008). This means that the turbulent velocity disper- 
sion is closely tied up with the gas surface density, which 
in turn is closely tied up with the star formation rate per 
unit area E^. In fact Genzel et al. (2011), compiled a large 
sample of the higher-redshift data described above to 
show that all galaxies with criD,Ha ^ 35 km/s have high 
star formation rate densities > O.IM© yr~^ kpc~^. 
At the same time, outflows are ubiquitous in all galaxies 
in which the star formation rate per unit area exceeds 
this same value of O.IM© yr~^ kpc~^ (Heckman 2002). 
Indeed, the star forming nucleus in M82, the archety- 
pal example of an outflowing starburst, has a star for- 
mation surface density of ~ 5Mq yr~^ kpc~^ (Strickland 
& Heckman 2009) and an average velocity dispersion of 

90 =b 30 km/s (Westmoquette et al. 2009). We may be 
just beginning to uncover the role of turbulent runaway 
in the history of galaxy formation. 
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